PROCEEDINGS OF THE 31 st ICRC, LODZ 2009 



1 



Atmospheric MUons from PArametric formulas: 
a fast GEnerator for neutrino telescopes (MUPAGE) 

M. Bazzotti*, S. BiagH, G. CarminatH, S. CecchinH, 
T. Chiarusit, A. Margiotta**, M. SiolH and M. Spurio*t 

* Dipartimento di Fisica dell'Universita di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy 
t INFN, Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy 
^INAF-IASF, Via Gobetti 101, 40129 Bologna, Italy 



Abstract. Neutrino telescopes are opening new op- 
portunities in observational high energy astrophysics. 
In these detectors, atmospheric muons from primary 
cosmic ray interactions in the atmosphere play an 
important role. They provide the most abundant 
source of events for calibration and for testing the 
reconstruction algorithms. On the other hand, they 
represent the major background source. 

The simulation of a statistically significant number 
of muons in large volume neutrino telescopes requires 
a big effort in terms of computing time. Some 
parameterizations are currently available, but they 
do not explicitly take into account the arrival of 
muons in bundles. A fast Monte Carlo generator 
(MUPAGE) was developed to generate single and 
multiple atmospheric muon events in underwater/ice 
neutrino telescopes. The code reduces the comput- 
ing time for the simulation of atmospheric muons 
significantly. 

The event kinematics is produced on the surface 
of a user-defined cylinder, virtually surrounding 
the detector volume. The flux of muon bundles 
at different depths and zenith angles, the lateral 
spread and the energy spectrum of the muons in the 
bundles are based on parametric formulas CD, 0. 
These formulas were obtained according to a specific 
primary cosmic ray flux model and constrained by 
the measurements of the muon flux in the MACRO 
underground experiment |3j|. 

Some MUPAGE applications are presented. 

Keywords: Simulation of atmospheric muons, neu- 
trino telescopes, Monte Carlo generator 

I. Introduction 

The most abundant signals for a neutrino telescope are 
due to high energy muons resulting from the extensive 
air showers produced by interactions between primary 
cosmic rays (CRs) and atmospheric nuclei. Although 
neutrino telescopes are located at large depths under 
water or in ice, taking advantage of the shielding effect 
offered by these detector media, a large flux of high 
energy atmospheric muons can reach the active volume 
of the detectors. 

The atmospheric muons represent an insidious back- 
ground for track reconstruction as their Cherenkov light 



can mimic fake upward going tracks. This kind of 
signatures can be confused with the cosmic neutrino 
signals that neutrino telescopes are searching for. On 
the other hand, atmospheric muons are a useful tool to 
test offline analysis software, to check the understanding 
of the detector and to estimate systematic uncertainties. 
Moreover, the pointing capability of the telescopes can 
be studied using atmospheric muons through the mea- 
surement of the moon shadow. 

Atmospheric muon bundles can be accurately repro- 
duced by a full Monte Carlo (MC) simulation (e.g. 
CORSIKA [4 1), starting from primary CR interactions 
with atmospheric nuclei, generating the resulting large 
number of air showers and propagating the muons 
through sea water or ice. A full MC requires a large 
amount of computing time and therefore the production 
of a large statistical event sample cannot be easily 
obtained for a detector of a cubic kilometer size. In a 
recent search for a diffuse flux of neutrinos [5|, the use 
of a full MC simulation limited the number of generated 
atmospheric muons at 63 equivalent detector days to 
represent a background for 807 active days of data. 

In order to save the big effort of computing time, a 
fast Monte Carlo generator is therefore essential. Some 
parameterizations for the underwater/ice flux and energy 
spectrum are available in literature 0, Q, (8). None of 
them, however, gives the possibility to simulate two or 
more muons produced in the same CR interaction (muon 
bundle). In this paper the event generator MUPAGE is 
presented. It is based on parametric formulas and takes 
into account the multiplicity of muon bundles. 

II. Simulation of atmospheric muon bundles 

Parametric formulas are derived in (TJ, |2| and de- 
scribe the flux, the angular distribution and the energy 
spectrum of underwater/ice muon bundles. The muon 
flux and energy spectrum are parameterized, for the first 
time, in terms of the bundle multiplicity m. 

From these parametric formulas an event generator 
called MUPAGE (MUon GEnerator from PArametric 
formulas) [9j was developed. The program output is 
a formatted ASCII table containing the kinematics of 
atmospheric muon bundles on the surface of a can. 
The can (Fig. [T| is an imaginary cylinder surrounding 
the active volume of any proposed detector between 
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1.5 and 5.0 km of water equivalent depth. The output 
table with the event kinematics can be used as input 
for the following steps of the detector-dependent MC 
simulation, which include the production of Cherenkov 
light in water or ice and the simulation of the signal in 
the detection devices as photomultiplier tubes. 




Fig. 1: Sketch of some input parameters. The cylinder sur- 
rounding the instrumental volume is the can, with radius J? C an 
and height i/ C an- The events are generated on an extended can 
with iicxt- The origin of the coordinate system does not have 
to be located at the center of the detector. The lower disk is 
at a depth -ff max with respect to the sea/ice surface. 

The code is available in a tar archive [9 1 or by sending 
a request to the authors. 

III. MUPAGE STRUCTURE 

The MUPAGE code is written in C++ and it has been 
tested with gcc version 3.2.x., 3.4.x and 4.1.x. The pro- 
gram requires ROOT libraries ffOl for the pseudorandom 
number generator. The tar archive contains the code, the 
Makefile, a README file, two template script files 
(for tcsh and for bash) to launch the executable and 
some subdirectories. After compilation, the executable 
and a Linux/ directory, containing the object files, are 
created. The core classes are contained in directories 
inc/ and src/, which do not need to be changed by 
the user. 

The simplest way to execute MUPAGE is to use, for 
tcsh, the C-shell script file run-mupage . csh, or, for 
bash, the Shell script file run-mupage. sh. Here, the 
user can modify the random seed, the run number and 
the number of events to be generated N gen . 

The flowchart of the Monte Carlo program is shown 
in Fig. [2] The bundle multiplicity, direction and impact 
point of the shower axis on the can surface are gener- 
ated first. According to the depth-zenith angle-intensity 
relation, this event can be accepted or rejected, using a 
standard Hit-or-Miss method. Once selected, if the event 
is a single muon, the energy is sampled from the energy 



distribution expected for that specific depth and zenith 
angle. If the event is a multiple muon, for each muon 
the radial distance R from the shower axis is sampled, 
followed by the muon energy (which depends also on 
R). The impact point of each muon in the bundle is 
the projection on the can surface of its position on the 
plane perpendicular to the shower axis, where it has been 
generated. Consequently, since all muons in the bundle 
are assumed to reach that plane at the same time, also 
the arrival time of each muon on the can surface is 
computed. Details are given in |9|. Since the multiplicity 
to of the bundle is generated according to the muon 
flux, it can happen that some muons in the bundle do 
not geometrically intercept the can surface and therefore 
they are not written in the output file. 
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Fig. 2: The flowchart of the MUPAGE event generator. The 
smooth-angle rectangles indicate the extraction of uniformly 
distributed random values. The decisional rhombuses in bold 
select values according to formulas reported in 1 1 1, with a Hit- 
or-Miss method. The procedure is iterated for N gen events. 

A. Description of the input file 

The generator needs some parameters as input: the 
dimensions and depth of the can surrounding the de- 
tector (see Fig. [TJ, the density of the detector medium 
and the ranges of the various simulation parameters (e.g. 
multiplicity, zenith angle, muon energy). The user can 
change the default values (which refer to a production 
done for the ANTARES experiment ifTTl ). The main 
parameters are: 

Hmax (2.475 km): vertical height of the sea/ice level 
with respect to the can lower disk; 
Zmin (—278.15 m): position of the can lower disk in 
the detector frame, Zmax (313.97 m): position of the 
can upper disk. The can height is defined as: H can — 
7 — 7, . ■ 
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Rcan (238.61 m): can radius; EnlargedCanR (100 m): 
can radius extension; 

density (1.03 g cm -3 ): mean value of the density of 
the detector medium (ice or water) to convert Hmax in 
km.w.e.; 

THETAmin (0°), THETAmax (85°): minimum and 
maximum value of the zenith angle range; 
MULTmin (1), MULTmax (100): minimum and max- 
imum value of bundle multiplicity range; 
Rmax (100 m): maximum muon lateral spread with 
respect to the shower axis; 

Ethreshold (0.001 TeV): the minimum energy threshold 
of all muons combined in the event bundle. 

The radius of the generation volume is increased by 
the quantity EnlargedCanR in order to accept also 
muons very far from the shower axis. Since one of the 
aims of the code is to save computing time, the user 
is encouraged to set the same value for EnlargedCanR 
and Rmax. 

An example of usage of parameter Ethreshold is 



given in Sect. IV-A 



B. Description of the output files 

The output file ($OUTFILEl in the template script) 
contains all information about the generated events in a 
formatted ASCII table. An example of an output data 
file with 1000 generated events (run_01 . evt) can be 
found in the tar archive. Each line of the table contains 
the following information: 

evt_id m track_id Xi yi Zi v x v y v z Ei ti G_id 
where: 

• evt_id is the event number; 

• m is the multiplicity of the muon bundle at the 
depth where the shower axis hits the can; 

% for each muon m c in the bundle intercepting the 
can (m c < m): 

- track_id — i (i = 1 , m c ) is the muon identifier 
in the event; 

- (xi,yi,Zi) are the coordinates [in m] of the 
muon impact point on the can surface; 

- (v x ,v y ,v z ) are the direction cosines of the 
muon, coincident with those of the bundle axis; 

- Ei is the energy [in GeV] of the muon; 

- ti is the time delay [in ns] of the i-th muon at 
the can surface with respect to the first muon 
(?' = 1). ti can be either positive or negative; 

- G_id is the GEANT particle identification 
number (6 = \i~ by default). 

The livetime and its statistical error (see Sect. |HI-C| > 
corresponding to the total number of simulated events 



N 



gen 



is written in a second file ($OUTFILE2). 



Note that there is no bias or ordering (in energy, 
multiplicity and zenith angle) in the simulation. All 
events have the same weight so the output file reproduces 
a real data file. A disadvantage of this approach is that 
since all the input physics parameters (primary CR flux, 
interaction models, muon transport) are fixed, the events 
cannot be re-weighted. 



C. Equivalent livetime of the simulation 

From N gen , the equivalent livetime of the run is 
automatically calculated. Since the flux of muon bundles 
of a given multiplicity m is known for each value of 
depth H at a given zenith angle di, the expected rate on 
the can upper disk at the depth H = H m [ n is: 

N MC (An t ) = t>(H min ,$ tl m) -S-AQi [s" 1 ] 

where, Af2j = 27r(cos$ii — cosi?2i) is the solid angle 
centered at $j = (flu + •&2i)/% and S — TrRl xt cos'&i 
[m 2 ] is the projected area of the upper disk. The equiva- 
lent livetime is computed from the number of generated 
events N(A£li) on the upper disk in the solid angle Af2j 
as: 

T(Afii) = N(ASli)/N M c(Mk) [s] 

The livetime is computed as the weighted average on 
33 different solid angle regions of T(Af2j), which have 
the same value, within statistical errors, and are written 
in the $OUTFILE2 file. 

IV. Examples and application 

The development of MUFAGE was motivated by the 
need of a large sample of atmospheric muons in order 
to simulate the response and the possible background 
for the neutrino studies in the ANTARES and NEMO 
Experiments and moreover to study possible detector 
geometries for the future cubic kilometer detector in 
the Mediterranean Sea, on behalf of the KM3NeT Con- 
sortium. Although a full Monte Carlo simulation is 
also available fl2l . MUFAGE offers the possibility of 
a faster simulation. Some examples of the MUFAGE 
performance are presented. 

A. Usage in the ANTARES Collaboration 

The ANTARES neutrino telescope ifTTI has been 
taking data since March 2006 and it has been com- 
pleted in May 2008. During the first 6 months of 
operation, in its one line configuration, measurements 
of atmospheric muons from data are taken. For this 
first analysis, two samples of atmospheric muon bundles 
with different multiplicity ranges are generated using 
MUPAGE: for to = 1 - 30, the 8.8 x 10 s generated 
events correspond to a detector livetime of 11.7 days; 
and for m = 31 — 1000, the 8 x 10 6 generated events 
correspond to 32.4 equivalent detector days. After the 
MUPAGE kinematics generation, muons are propagated 
inside the can, where the emission of Cherenkov light 
and the production of signals in the photomultiplier 
tubes (PMTs) are simulated with a software package lfl2l 
which needs a computing time almost a factor 10 larger 
than that required by MUPAGE. The two samples are 
used to investigate the efficiency of the track recon- 
struction method, reproducing the time and the angular 
residuals very well, and to measure the vertical muon 
intensity versus depth [ 13 1, in good agreement with other 
published values. In particular for the latter goal, the 
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MUPAGE simulation was used to convert the measured 
rate of reconstructed tracks to the single muon intensity. 

For the 5 line configuration of the ANTARES detector, 
in 2007, MUPAGE is used to generate a data set with 
livetime of one month. This production is split into 322 
files (< 2 GB each). Each file contains 10 7 events, 
corresponding to 8060±3 s. The computing time needed 
to produce a file (on a 2xlntel Xeon Quad Core, 2.33 
GHz per core) is less than 2 hours. This MUPAGE data 
set is used to compare the zenith and azimuth angle 
distributions of reconstructed tracks with the real data 
and with a CORSIKA full Monte Carlo simulation fl4l . 
A preliminary comparison between data and MC of 
electromagnetic showers produced by high momentum 
muons is also performed [15|. The MUPAGE simulation 
enables the development of a deconvolution method 
which to be applied to the real data in order to obtain 
the experimental atmospheric muon flux as a function 
of the sea depth |[T6l . 

The same MC data set is additionally used to estimate 
the background rate induced by simultaneous muon 
bundles originating from different cosmic rays. If two 
atmospheric muons arrive during a trigger window, the 
produced signals on the PMTs cannot be distinguished 
from each other. In particular, the timing patterns of 
the light from the two tracks can be such that the 
reconstruction result is a single upward going track. 
The rate of triggering 'coupled events', which reach 
the 5 line ANTARES detector within a time window 
smaller than 4 /its, is about 7 x 10 -4 Hz. However, no 
events are reconstructed as upward going, fixing the 90% 
confidence level upper limit at 9 x 10~ 7 Hz. 

A third data set is produced for the background study 
of the high energy neutrino diffuse flux (above 100 TeV). 
In order to optimize the computing time, a (conservative) 
cut on the minimum energy threshold of all muons 
combined in the bundle is applied, E t h r > 3 TeV, and 
the event multiplicity is divided in several ranges, split 
is 852 files. The total computing time is 232 hours with 
a livetime of one year. 

B. Usage in the NEMO Collaboration 

The NEMO Collaboration El has looked for an op- 
timal site for the installation as well as for the develop- 
ment of technologies for an underwater cubic kilometer 
experiment. The NEMO Phase-1 project, operating from 
December 2006 to May 2007, has allowed a first valida- 
tion of the feasibility of cubic kilometer detector (TH). 
A small MUPAGE MC data set is produced, generating 
4x 10 7 events corresponding to ~ 11.3 equivalent hours. 
The MC simulation of zenith and azimuth distribution of 
reconstructed tracks is used for the comparison with 1 1.3 
hours of data taking, giving an excellent agreement |fl9l . 

C. Usage in the KM3NeT Consortium 

The KM3NeT Consortium aims at the definition of 
a complete project for an underwater cubic kilometer 
neutrino telescope to be installed in the Mediterranean 



Sea ll20l . Since the final detector performance will 
depend critically on the rejection of atmospheric muons, 
MC simulations for different detector geometries as a 
function of site depth are performed. For example, a 
small MUPAGE MC data set with a livetime of 105 min 
is used to study the response to atmospheric muons for 
different PMT orientations EH 

V. Conclusion 

MUPAGE is a fast generator of the kinematics of 
atmospheric muon bundles. It uses parametric formulas, 
valid in the interval 1.5 < h < 5.0 km water equivalent 
and ■& < 85°, for the flux of single and, for the first time, 
multiple muons. The simulation of the energy spectrum 
of single and multiple muons takes into account the de- 
pendence of the muon energy on the shower multiplicity 
and the distance of the muon from the shower axis. 
The generator, developed to minimize the computing 
time, represents a useful tool for underwater/ice neutrino 
telescopes for the production of a large amount of 
simulated data. The code is currently used for the Monte 
Carlo simulation of atmospheric muons on behalf of the 
ANTARES and NEMO Collaborations. Furthermore, the 
KM3NeT Consortium is using the code to study the 
atmospheric muon rejection performances of different 
detector geometries. 
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